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A FAST-MARCHING ALGORITHM FOR NON-MONOTONIGALLY 

EVOLVING FRONTS * 

ALEXANDRA TCHENG, JEAN-CHRISTOPHE NAVE t 

Abstract. The non-monotonic propagation of fronts is considered. When the speed function 
F : MF X [0,T] ^ M is prescribed, the non-linear advection equation (j)t + TlV^I = 0 is a Hamilton- 
Jacobi equation known as the level-set equation. It is argued that a small enough neighbourhood of 
the zero-level-set Ai of the solution cj) : x [0, T] ^ M is the graph of ijj : M where solves a 

Dirichlet problem of the form if(u, V^(u), V'0(u)) = 0. A fast-marching algorithm is presented where 
each point is computed using a discretization of such a Dirichlet problem, with no restrictions on 
the sign of F. The output is a directed graph whose vertices evenly sample At. The convergence, 
consistency and stability of the scheme are addressed. Bounds on the computational complexity are 
estimated, and experimentally shown to be on par with the Fast Marching Method. Examples are 
presented where the algorithm is shown to be globally first-order accurate. The complexities and 
accuracies observed are independent of the monotonicity of the evolution. 

Key words, front propagation, Hamilton-Jacobi equations, viscosity solutions, fast marching 
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1. Introduction. Front propagation is a time-dependent phenomenon occurring 
when the boundary between two distinct regions of space is evolving. It is possible to 
make the distinction between monotone and non-monotone motion of fronts. Consider 
a fire propagating through a forest: the interface divides space into a burnt and an 
unburnt region. It evolves monotonically in that, if a point x = (xi, ... ,Xn) of space 
belongs to the burnt region, then it cannot belong to the unburnt region at a later 
time 20 . In contrast, if a plate containing water is gently rocked, the front separating 


the dry region from the wet one may advance or recede. In this paper, we present an 
algorithm that dynamically builds a sampling of the subset of x [0, T) consisting 
of the surface traced out by the front as it evolves through time. The structure of 
the scheme is akin to the Fast Marching Method, yet our approach is oblivious to the 
monotonicity of the evolution. 

Given an initial front Co as a codimension-one subset of and an advection rule, 
our goal is to recover the front Ct at later times t > 0. Let each point of the front 
evolve with a prescribed speed F : x [0, T] ^ M in the direction of the outward 

normal to the front vit : Ct ^ S^. The evolution is non-linear, such that even if F 
and Co are smooth the front may become and undergo topological changes |^. 

On the one hand, a robust numerical method for tracking either kind of evolution 
is the Level-Set Method (LSM) 22,^. This implicit approach embeds the front 
as the zero-level-set of a function 0 : x [0, T) ^ M, and solves the initial-value 

problem: 


r + = 0 inM-x(0,T) 

\ 0(x, 0) = 00 (x) on X {0} ^ ^ 

where 0o satisfies {x : 0o(x) = 0} = Co- Assuming the computational grid comprises 
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points, the total complexity of the first-order algorithm solving is 0(Ar^+i) 
when endowed with the usual CFL condition At oc 1/N. On the other hand, the Fast 
Marching Method (FMM) 1^-28 33 may be used when F = F{x.) > S > 0 and is 


therefore suited for monotone propagation. This approach reduces the dimensionality 
of the problem by building the first arrival time function -0 : ^ M, i.e., t = 

gives the unique time at which the front reaches x. This function solves the boundary- 
value problem: 


ly^l = ^ in So 

'^(x) = 0 on Co 


( 1 . 2 ) 


where So is the unbounded component of \Co for F > 0. The FMM solves (1.2) in 
0{N^ \ogN^) time using a variant of Dijkstra’s algorithm |^. Schemes derived from 
the FMM include: an extension to the case where F depends on time [M,35 ; and 
the Generalized FMM [^. Problem (1.2) may also be solved using the Fast Sweeping 


Method 1^, which is a first-order accurate method running in 0{N^) time. When 
F = 1, higher order methods exist, such as the one of Cheng and Tsai which 
exploits the relation between and the time-dependent Eikonal equation. 

The GFMM of allows F to vanish and change sign albeit at the expense of 
the computational time. In [^, when n = 2 the authors propose a scheme that can 
handle such speed functions while featuring a computational complexity comparable 
to that of the FMM. This approach considers the set M := {(x,t) : x G Ct}. Each 
point p e M may be described as p = (x, p,'0(x, p)) and/or p = ('0(p, t), p, t) and/or 
p = (x,'0(x, t), t). In regions where |F| > > 0, the algorithm builds using the 

EMM, whereas near points where F = 0, it solves a PDE satisfied by either or fj. 
Despite a number of satisfactory results, the mechanism to change representation is 
cumbersome, and A4 is undersampled near regions where F vanishes. 

In the present paper, rather than limiting ourselves to the n + 1 representations 
just mentioned, we now describe A4 locally with functions of the form '^(u) where 
the points u = (ui, ... ,Un) belong to some hyperplane lying in x [0,T]. Consider 
the orthonormal spanning set for this hyperplane consisting of {lii ,... ^Un}^ and let 
the normal be fin+i- So that in this convention, any point (u, r^n+i) of space-time 
MF X [0,T] may be written as (u,= uiui + ... + UnUn + Section]^ 

shows that M solves a Dirichlet problem of the form: 


F(u,'0(u), V'0(u)) = 0 in Z// C 

'0(u(x,t)) = ?Xn+i(x,t) on xGCtHV 


(1.3) 


for appropriate neighbourhoods U and V. The Hamiltonian function H : x M x 

^ M depends on through the speed function F and involves constants that 
capture the relative orientations of the two coordinate systems in use. Section 
presents the discretization of (1.3) using finite-differences, as well as the design of 
a constrained minimization problem. Section discusses the different parts of the 
algorithm for the case n = 2. The protocol we propose is similar to the EMM: Points 
sampling M first belong to the narrow band JV before moving to the accepted set A. 
When Pa G A/" is accepted, a local {ui ,..., lin+il-coordinate system is found. Using pa 
and another point in M, a new point is computed through the finite-difference solvers 
and the optimization problem. Section highlights the properties of the solvers and 
the global scheme. In particular, convergence of the local solvers is addressed and the 
total computational complexity is estimated. As is illustrated in ^with numerous 
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examples, the output of the algorithm is a directed graph. Its meshsize is bounded 
below by a predetermined parameter h, and its vertices provide a discrete sampling 
of Ad. Given t G ( 0 ,T), the front Ct can be recovered using interpolation. 

Methods with a similar flavour have been explored in 
order fast interface tracking method, and 


17 


24 


which presents a high 
where the authors formulate n ‘quasi- 
linear PDEs satisfied by the manifold’s local parametrization’ and ‘solve that system 
locally in an Eulerian framework’. 

Eor clarity of illustration, the setting of our paper is n = 2, so that x = (xi, X 2 ) = 
(x,^), u = {ui^U 2 ) = ('U, u) and us = w. It is worth noting however that most of the 
underlying ideas, which are presented in ^and ^ easily extend to the general higher 
dimensional case. Current obstacles to the full generalisation of the method lie in the 
design of a practical interface tracker, able to capture the global features of the front 
at all times. The successful results of Sections]^ andwhen n = 2 provide a proof 
of concept that such a goal should be pursued in future work. 

Our approach offers numerous advantages when n = 2. It extends previous work: 
If the local and global systems coincide then -0 = -0, whereas if {ui^U 2 ) = (^, t) and 
Us = X then ip = ip. By construction, the transition from one representation to another 
is smooth, and the resolution of the sampling is regular. The initialization is almost 
identical to the main procedure, and does not require any information away from Cq. 
The accuracy of the scheme is 0{h). If Cq is sampled by m points, the computational 
complexity is bounded by CN where C = max{0(m), O (l0’^+^)} with N = 0{mA). 
When run on a monotone example, for the same computational time, our method 
yields a more accurate solution than the EMM. 


2. Hyperplane representation. In this section, we derive Problem (1.3) from 
the IVP O) after making a few assumptions. We then argue that the solution 
of ( |1.3| ) locally describes Al, and that a finite number of such solutions provides a 
covering of Ad. 


2.1. Assumptions. The set Co is a closed, co-dimension-one subset of C 
without boundary, where O is bounded. Moreover, Cq is the graph of a function. 
The speed function F : x [0,T] ^ M is continuous. It may vanish and change 

sign. Under those assumptions, the LSE in is a Hamilton-Jacobi equation with 
a unique continuous viscosity solution [ 6 ||^. It follows that Ad embeds in x [ 0 ,T] 
as a manifold. We denote as z^(x, t) the outward normal to Ad at (x, t). 


2.2. PDE on an arbitrary plane. Let a normal vector u E he given, and 
consider a corresponding plane lying in xyt-space. We denote as uvw the unique right- 
handed orthonormal coordinate system satisfying the following condit ions : 4) = P, the 
hu)-plane is vertical and the t-component of v is positive. See Eigure 


2.1 


We have: 


U } 

\ 1 

( CXl 
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V 
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w j 

f ) 

V 71 

72 

73 



( 2 . 1 ) 


where the o’s, /3’s and 7 ’s are constants. In particular, 03 = 0 and Ps > h (see 
Appendix [A| for details). Let us now assume that a neighbourhood V of (xo,to) G Ad 
may be represented as the graph of the function: 

^ M p; :u^ p;{u) ( 2 . 2 ) 


If u • z>(xo, to) > 0, then by the Implicit Eunction Theorem 0(x, t) = w — '^(u) where 
(p is the solution of Using implicit differentiation, we replace the derivatives 
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V 



(b) Two neighbourhoods and their 

(a) A4 C R X [0,T]. The xj/i-coordinate system is ^^^.coordinate systems, 
the global one. 


Fig. 2.1: Illustration of the notation for n = 2. 


appearing in Oby 

(f>xi = 7 i - «iV’« - Mv « = 1, 2 and (/>* = 73 - aaV’w - /?3V’« (2.3) 


Using the orthonormality of the ixure-coordinate system, (1.1) can be rewritten solely 
in terms of 0 ^ 3 , /Ss and 73 as 

[7 - (aipu + Pi^v)] + G\J+ - {aipu + I3ipv)f = 0 (2.4) 

where the subscript -3 was dropped, and G(u,'^(u)) := F (x(u,'0(u)), t(u,'0(u))) 
was introduced. Equivalently, Equation \2A\ can be rearranged using the vectors 
R := {a, P, 7 ) and u := {-ipu, 1), as: 




^u • u — (^ly ' R^ 


= 0 


(2.5) 


Note that this formulation is independent of dimension. In the event where u = 
^(xq, to), the point (xq, to) is a local maximum of V, and we may explicitly relate the 
orientation coefficients to the speed function: 


F= 0, 


-F(xo,to) 


( 2 . 6 ) 


a/1 + F2(xo,io) ’ x/l + F2(xo,to )J 
We arrived at the following conclusion: Suppose 7 /^ satisfies the boundary condition: 

'0(u(x, t^)) = u;(x, t^) when x G C^o fl V (2.7) 


and solves (2.5) on U C M?. Then letting t* := t(u,'0(u)) for u G tt, we have that 
X* := x(u,'0(u)) G Ct-. 
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2.3. Well-posedness. Equation (2.4) is a Hamilton-Jacobi equation of the form: 


:= -[7-/3V^^]-G(u,V^)yV'2 =0 (2.8) 

where we defined the Hamiltonian x M x ^ M. The well-posedness of this 

equation is guaranteed in the class of viscosity solutions [T . The following theorem 
justifies defining H as the negative of the left-hand side of ( 2.4| ). 

Theorem 2.1. The viscosity solutions of Equations and are equivalent 
in the sense that the zero level-set of (j) at various times t = t(u, w) is precisely the 
set of points x = x(u, vu) for which '^(u) = w. 

Note that this equivalence is investigated in the general cas e in O sher’s paper 21 . 
Proof If A4 is at a point p, then the manipulations of ^ 2.2 hold in the strong 


sense. Thus, let j\4 be singular at the point p, and consider the second-order PDE: 


=eA(/)^ e >0 


(2.9) 


i.e., A viscous term was added to the right-hand side of the level-set equation. Define 
:= w — (j)^{u^v). Eollowing the same reasoning as in ^2.2, using the orthonormal 
relations, as well as the fact that 0^3 = 0, Equation (2.9) becomes: 


[7 - + G\J 7/^2 + (^ + + {Pi + Pi) i^lv] (2.10) 

Taking e to zero, the arguments detailed for example in §1.5 of Lions’ book apply 
to yield that the Hamilton-Jacobi equation of interest is ( |2.8[ ). □ 

In the light of this result. Equation (2.8) is preferred over (2.4) in the rest of this 
paper. 

2.4. Geometric properties. We illustrate how M can be globally described 
by a covering of T/^-functions. 


ip provides a local representation of M. Suppose that the solution of is such 
that V \ = U xfj{U) is G^. Let tmin = inf{t : (x, t) G V} and t^ax = sup{t : (x, t) G V}. 
Theorem 2.2. For any given t* G (tmin, ^max); we have: 


{x G : '0(u(x, t*)) = u;(x, t*)} = Ct* fl V 


The proof follows the same argument as the proof of Theorem 4.2 in 
therefore omit it. 


30 , and we 


Covering of M. It follows from compactness that M can be covered by a finite 
number of images of functions -0, each solving a problem of the form (1.3). The 
construction presented below builds one such covering that is particularly suited for 
algorithmic purposes. Given At > 0 and 0 < to < T — At, consider the strip 


St := {(x, t) G A4 : to < t < to + At} 


( 2 . 11 ) 


Theorem 2.3. Suppose thatCt^ has codimension one inW^ x [0,T). There exists 
a finite covering of St consisting of the graphs of functions ipi, where each ipi solves a 
Dirichlet problem of the form ( fi.7| ) with boundary conditions imposed on Ct ^. 

Proof Consider an open covering of Ciq by open segments s with the following 
property: Defining 

= — / t)ds where \s\ = I ds 




-t 


( 2 . 12 ) 
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for each p E s we have: i>{p) • z/q^ > 0 . Since Ct^ is compact, we may extract a finite 
open covering, say Sq = : i = 1,2 ,...,/}. Fixi, and define the future domain of 

influenee of Si as the subset of St satisfying: 


Wi :={(x, t) : 3{(xn,tn)}^i ^ (x, t)such that for each 

(x^,t^) 3 a characteristic of the LSE starting in Si and ending at (x^,^^)} 

(See for example for details on the method of characteristics.) See Figure [2^ for 
an illustration. We have St = Define: 

A't := J C'dWi where \Wi \ = J dWi (2.13) 


Note that limAt^o Picking At > 0 small enough, it follows from the 

continuity of A4 that for each p G Wi, we have: z>(p) • > 0. Consider the Dirichlet 

problem i/(u, VV^;/t(z/^^)) = 0 with boundary conditions imposed at Si C Ciq. 
From Theorem |2.2[ Wi is contained in the graph oi fj. U 

A global covering may then be obtained from decomposing M into strips of the 
form (2.11) and extracting a finite subcover. 


3. Discretization. Throughout this section, assume that the iz'zzrc-coordinate 
system is fixed. Suppose we are given two points belonging to M of the form 
Pa = (ua,'^na) and pi) = (u5,rc5); and wish to compute Wd = f^{ud) at a location 
Ud = {ud^ '^d)- We will refer to Pa and pi^ as the parents of the ehild point pd = (ud, Wd). 
In §3.1[ we first assume that Ud is given and propose two different solvers that take 
as input Pa^ Pb and Ud, and return a value Wd- Then in §3.2[ we turn to the ques¬ 
tion of how Ud should be determined. The answer takes the form of a eonstrained 
minimization problem^ for which we propose two different methods. 

Section thus provides a framework for solving the problem locally. Global issues 
also need to be addressed - See §4.4[ 

3.1. Solving the PDE. As illustrated in Figure [tT| define: 


Sa := Ud - Ua Sb := Ud - Ub Sc := Ud - Uc 


(3.1) 


and let Si = ((si)u, (si)v) where i = a, b, c he the corresponding unit vectors. Then 


'ipi = 


Wd - Wi 


where i = a^b^c 


(3.2) 
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Fig. 3.1: The points Pa = {ua,Va,Wa), Pb = {ub,Vb,Wb) and Pc = {uc,Vc,Wc) belong 
to A. Given = (ud^Vd)^ we wish to compute Wd- 


provide order approximations of the directional derivatives of ip at Making 
use of the calculus identity: for any unit vector k, we have: 


f 

V 


(<^a)u (^Sa)v A f V^a A ^ — 1 ( 

('^i)u {^^i)v J \ J \ 


where i = b,c (3.3) 


provided that Sa and Si are not colinear. We write 

i> = —Mwd — N with M := (m^i, 0) and TV := (n^^, —1) (3.4) 


where the constants rriu^ riu and Uy depend on (ua^Wa), (u^,re^) where i = b or 
c, and Ud. The quantity D provides a order approximation of u. 

3.1.1. Direct solver. Consider approximating G by setting it equal to Go = 
G{pa)- Note that this is independent of Wd, which implies that the relation 


u ' R Gq 




= 0 


can be rearranged as a quadratic in Wd^ Indeed, first rewriting (3.5) as 

2 


^ (l + Gq) — Gq (z> • k) 


(3.5) 


(3.6) 


and then making use of (3.4), we arrive at 
, 2 


(^i> • = kiw^ + ‘^k 2 Wd + ks and z> • i> = + ^6 (3-7) 


where 


ki = (^ R-My k 2 = (^R ■ (^R ■ fc3=(E-7v)^ 

ki = M -M k 5 =M -N ke = N-N 

Rearranging further yields 

(/ci + GQ{ki — k^)) Wd-\-2 (/c2 + G‘^{k2 — k^)) Wd + (fe + Go(A:3 “ ^e)) = 0 


(3.8) 
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The discriminant of this quadratic reads: 


0 + PiGq + P2Gq — Gq (pi + P2Gq) (3-9) 

where pi = —2k2k^ + kikQ + /C 3 /C 4 and P 2 = Pi + ^5 ~ k^kQ. We set 

- [k2 + Go(fe “ ^5)) + Go Vpi + P2GI 

-71- , ^2/7 -TTl- 

(/ci + Gq(/ci — /C 4 )) 


Since sign( 7 ) = —sign(Go), this choice of root minimizes td = aud + ^Vd -\-^Wd where 
a = 0 and /3 > 0. This is consistent with the control theoretical approach to this 
problem [TTj[l^[33|[34| . If 0 = /ci + G^ki — /C 4 ), the quadratic relation degenerates to 
a linear one with solution: 


ks^Gliks-ke) 

2ik2^Gl{k2-k,)) 


(3.11) 


The outward normal to M. at pd can then be c omputed as Pd = k/\P\. 

Remark: As was already mentioned in ^ 2.3, when R = (0, 0, ±1), the PDE reduces 
to the Eikonal equation. If working on a Cartesian grid, formulas (3.10) and (3.11) 
can be verified to agree with the solvers used in the traditional EMM. 

3.1.2. Iterative solver. An alternative to approximating the speed term G as 

(3.12) 


in ^3.1.1 is to solve (2.5) iteratively in pseudo-time r. i.e., Let 

v'^-.= -Mw2-N and G" := F( x(ud, Wd), i(ud, ) 


and iterate: 




P^ ' P^ — (^P'^ • 


w 


n+l 




At 


(3.13) 


until — w2\ is below some pre-determined tolerance, 

be used to initialize w^. 


Definition 5 of 19 


see 


The direct solver can 
The size of the pseudo-time step is determined based on 
Appendix [b] for details. 


3.2. Optimal sampling. Suppose that a set A of points sampling a part of At 
is available. Given Pa^ Pi G M, we now consider where to compute a new point pd- 
Constraints are heuristically imposed so as to efficiently build a regular sampling of 
At. Rigorous justifications of these choices are provided in ^ where the properties of 
the scheme are analyzed. 


3.2.1. Constrained minimization problem. 

Aeeuraey. According to (3.2), the accuracy of either solver increases as \sa\ and 
I Si I become smaller. We combine those requirements into a single function to be 
minimized: 


/(Ud) := |sap + |sy 

The result is that the child point pd should he close to its parents. However, we now 
argue that it should not be too close. 

Evenness of the sampling. Repulsion between pd and its parents is introduced 
to avoid oversampling parts of Al. The minimum three-dimensional (or (n + 1 )- 
dimensional) distance between pd and any other point in A should be at least h, 
where h is a predetermined parameter. As will be clear in ^ h bounds the meshsize 
of the resulting graph from below. 
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Characteristic structure. Lastly, we make sure that pd does not violate the char¬ 
acteristic structure of the solution: Wd must be real, and in addition pd has to satisfy 
td > max{ta,ti} for causality to hold. We impose the more stringent condition: 


td > max{ta, ti} ^ At where At = 


h 




with j = a if max{ta,ti} = ta and j = i otherwise. This choice of At speeds up 


computations and simplifies our analysis, as will be evident in ^ 5.1.1 where At = Ph. 

Optimal sampling - Summary. Given a pair of points Pa^ Pi ^ A, we wish to 
minimize the objective function 


subject to the constraints: 



(A) 

gi{ud) = ^{wd) = 0 

(VI) 

52(ud) = td- ma.x{ta,ti} > At 

(V2) 

gsiud) = min{||pd -p||} > h 

peA 

(E) 


3.2.2. Solving the constrained minimization problem. Although the ob¬ 
jective function is simple, the constraints are non-linear functions of Unlike the 
other two, constraint (VI) does not require evaluating Wd- If using the direct solver, it 
amounts to checking the sign of the discriminant. Constraint (V2) requires conversion 
of the data to xt-coordinates. Constraint © can be very costly to verify if A is large. 
This is why it is preferable to work with a predetermined small subset of A. To be 
efficient, a scheme solving this problem should not require too many evaluations of 
the constraints. We briefly discuss two possible methods, and relegate the details of 
the procedure to Appendix [C] 

The grid method. This simple yet efficient approach is the one we use in practice. 
The problem is solved using grids sampling a neighbourhood of Pa and Pi in the uv- 
plane. The objective function / is computed at those points where all three constraints 
are satisfied. The location of the minimum of / is found, and a finer grid centered at 
this point is defined. The procedure is repeated until convergence e.g.. The change in 
the value of the minimum is below some pre-defined tolerance. 

The Lagrangian method. Following (§2.1, 2.2, 3.1), the augmented Lagrangian 
L is defined by adding the three constraints to / along with Lagrange multipliers pi 
and penalty coefficients q > 0, where i = 1,2,3. The following procedure is iterated 
until convergence e.g.. The change in the value of the minimum is below some pre¬ 
defined tolerance. First, for some values of the q and the unconstrained minimum 
of L is found. This can be done using BFGS along with line minimization 23 . Then 


the q’s are increased and the are updated. Although convergence is guaranteed, 
this method requires a large number of evaluations of the constraints, which makes it 
too slow for our purposes. 


4. Algorithms. We present the algorithms used to generate the results in ^ 

4.1. Initialization. Little information is required to initialize the algorithm. 
The input is a sampling of Cq consisting of m points and m — 1 undirected segments. 
The normal z> at each point and a Final Time must also be provided. 
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Store those pieces of information into a list, say C. In the following, we will 
assume that each row of C contains data pertaining to one point, e.g., if p = (x,t) = 
(xi, X 2 , . •., t) has normal z>(x, t), then C might initially look like: 



Xi 

X2 


t 

vi (x, t) 

Z> 2 (x,t) 


7’„+i(x,t) 

row 1 

row 2 

-.1625 

-.1705 

0.0080 

0.0099 


0.0342 X lO-i'^ 
0.0398 X 10-15 

-.1634 

-.1706 

0.301 

0.385 


-0.5668 

-0.5909 


The variables Mj\f and help monitor the number of points in M and respec¬ 
tively. See Algorithmic 


Algorithm 1 Initialization 

1: A ^— >C, i — TTl T 1 

2: J\f i — Mj\f i — Tfl T 1 

3 : (minp,,e£ ||p-g||)/2. 


4.2. Get a local representation. This routine is called when a point Pa is 
accepted to go from a global to a local representation. See Algorithm 


Algorithm 2 Get a local representation 


1 

2 

3 

4 

5 


Find the L closest neighbours of Pa among A and form the list S. Set A = S. 
Find u such that u • i>{pi) > 0 for each pi e S. The default value is P = z/^. 
Compute the izizre-coordinates of the points in S using relation (2.1). 

Find the point Pb ^ S lying closest to Pa with Ub — Ua > h/2. 

Return: A, the change of coordinates matrix M, and pb. 


4.3. Computing a new point. Once the minimum is found using the direct 
solver, the iterative solver is run to improve the accuracy of the solution See 
Algorithm 


Algorithm 3 Compute a new point 


At = h/ A/G2(Ua) + 1, Pd = 0, i>d = 0. 

Solve the constrained optimization problem defined by p^, pi^ G(ua), h and At 


to get Ud. See §3.2| and Appendix [C] for details. 

if a solution to the optimization problem was found then 

Compute Wd associated with using the iterative solver detailed in ^ 3.1.2 
Compute the normal Od at {ud^Wd)- 

Get pd and i>d in (x, t)-coord. using the change of coordinates matrix M. 

end if 

Return: [p^; i>d]- 


4.4. Global constraints. The global features of the manifold are monitored 
with the help of the structure Book which consists of a list of segments. For example, 
if Af has the form: 

•V = [pi;p2; • • • ;PMa/-]^ then Boofc = [p2 -Ps; Pii -P23; • • •; Pi 4 -( 4 . 1 ) 
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Segments are not directed, i.e., [pi —P 2 ] is the same as [p 2 —pi]> Initially, the collection 
of segments lies in the x^-plane and provides a piecewise linear approximation of Cq. 
We impose the following constraints on Book. 

• Multiplicity: Each point p e JV appears exactly twice in Book. 

• Loops: A segment cannot start and end at the same point. A segment can only 
appear once in Book. 

• Intersections: Segments cannot intersect. 

• Spikes: The acute angle formed by two segments sharing an endpoint cannot be less 
than « 0.27r. 

The situations are illustrated with simple examples on Figures |4.1| 


4.4 The 


framed subfigures are obtained after connecting hanging nodes which are nodes that 
appear only once in Book. This is done as follows: 

Stitching hanging nodes. Given a hanging node pi = (xi,ti), let p 2 = (x 2 ,t 2 ) be 
the hanging node that minimizes the x^-distance ||xi — X 2 ||. Then add the segment 
[Pi — P2] to Book. Repeat until there are no hanging nodes left. 

Those global constraints are imposed to JV in Algorithm Note that we omit 
the details of the updating procedure of Book^ which consists in relabelling nodes and 
possibly deleting some segments. See 


29 


for details. 


Remark: When a new point is computed by Algorithm the checks performed 
at lines 7 and 8 need only involve the segments attached to that new point, which 
requires an 0(1) number of operations. 


Algorithm 4 Book keeping of Af 
1: if Mjv" > 1 then 

2: if a new point was computed by Algorithm then 

3: Update Book. 

4: else 

5: Update and clean Book. 

6: end if 

7: Check intersections, and clean Book if a point was removed from Af. 

8: Check spikes, and clean Book if a point was removed from Af. 

9: end if 


10: procedure Clean Book. 

11: do 

12: Multi = 0; Loop = 0; Inter = 0; 

13: Check multiplicity, and set Multi to 1 if a point was removed from JV. 

14: Check loops, and set Loop to 1 if a point was removed from Af. 

15: Check intersections, and set Inter to 1 if a point was removed from Af. 

16: while Multi+Loop+Inter > 0 

17: end procedure 
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Fig. 4.2: Loops. Once the double-loop is removed, the nodes are either stitched (top 
situation) or removed from M (bottom situation). 




Fig. 4.4: Spikes 
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4.5. Main loop. A couple of remarks on Algorithm are in order. Lines 5-7 
are ignored during the first m iterations of the Main Loop: This can be considered as 
the end of the initialization procedure. The fact that the algorithm does not compute 
a new value when ta > Final Time implies that JV is eventually empty, at which 
point the algorithm naturally ends. If Algorithm fails, next-to-nearest neighbours 
are used (line 14). Since \A\ < L, this occurs at most L times. Here, we present the 
algorithms with i = b. The procedure is unchanged if i = c. 


Algorithm 5 Main Loop 
1 : counter = 1 
2 : while Af is not empty do 

3: Find p G Af with the smallest t value, call it Pa = (x^, ta). 

4: Remove [pa] T'a] from Af. ^ — 1 

5: if counter> m then 

6 : Add [pa] T^a] to A in the M^-th row. ^ + 1 

7: end if 

8 : if ta < Final time then 

9: Set A = A and get a local representation using Algorithm 

10 : Run Algorithmwith the pair Pa-Pb to get pd^ i>d‘ 

11 : if Algorithm was successful at finding a new point then 

12 : Add [pd] T^d] fo Af in the Myv'-th row. ^ + 1 

13: else 

14: Remove [pb] i^b] from M, find a new pb G M, and go back to line 10. 

15: end if 

16: end if 

17: Run Algorithm 1^ to impose the global constraints on Af. 

18: counter ^ counter +1 

19: end while 


5. Properties of the scheme. We first discuss local properties of the scheme by 
approximating the solution of the optimization problem, and asserting the convergence 
of the solvers. We then turn to global properties. 


5.1. Local properties. Let us assume that Pa and pi £ A are exact, i.e., Pa^ 
Pi G A4; and that \\pa — Pi\\ = 0{h). Moreover, consider that we are working on a 
neighbourhood where -0 is so that i> = i>a and Wi = Wa + o{h?). We investigate 
the properties of the algorithms used to compute pd. We show that the constraints 
of the optimization problem are such that the discretizations of Equation (2.5) are 
monotone. 


5.1.1. Optimization. The objective function can be rewritten as: 

f{ud) = 2 (||ud - Uminll^ - Ua • uf) where Ui„i„ = (5.1) 

i.e., / is a paraboloid which reaches its minimum at Umin- It follows from the convexity 
of / that the solution of our optimization problem either is Umin, or lies on the 
boundary of the feasibility set. In the limit where h ^ 0, we expect Wmin ~ {wa + 
Wi)l2. Recalling the relation t = /3v-\-jw, if 7 is small we may have t^iin ^ max{ta, 4 } 
which violates constraint (V2). As a result, in general Umin does not belong to the 
feasibility set. 
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Turning to constraint (VI) requires the study of the sign of pi + P 2 G 0 . It can be 
shown that 


Pi > ||M|p||V|pe for some e = e(^,M,V) >0 (5.2) 

Picking (5 > 0 such that e > (^^/(l + (^^) gives pi + P 2 G 0 > 0 when |Go| < S. When 
|Go| > > 0, using the approximation Wi = Wa yields the following estimate: 

Pi + GIp 2 ^ {p - 0^Gl) ml + ml = > 0 (5.3) 


Next, we momentarily focus on the solution to the following subproblem: Mini¬ 
mize (A) subject only to constraint (V 2 ). The constraint is equivalent to 

P {Vd - fa) + 7 {Wd - Wa) > => {Vd - fa) + o(/l^) > h (5.4) 

SO that it primarily constrains Thus, as h ^ 0, the solution tends to 

Ud=Umin Urf = max{h + Ua,Umin} (5.5) 

The interplay between (V2) and (E) must be addressed. If pj G A since ta > tj, 

P {vd - Vj) + 7 (rcrf - Wj) > ph {vd - vj) + 0{h^) > h (5.6) 


and so \\pd — Pi\\ > h as h —> 0. 

In conclusion: Under the assumption that -0 is locally G^, in the limit where 
h —0, we have that (VI) holds, and satisfying constraint (V 2 ) implies satisfying 
constraint (E). It follows that the solution of the optimization problem is given by 
the minimum of subject to the constraint ( |5.4| ). Note that in the event where 
the assumption on the regularity of is dropped, neither existence nor uniqueness 
can easily be asserted. 


5.1.2. Direct solver. Given a location the direct solver returns a value 


dir 


using either (3.10) or (3.11). We first proceed to showing that this discretization is 
degenerate elliptic in the sense of [^, i.e.. Letting 


^f('*/’i,V’a) := -[7-/?V’«(V’a,'*/’i)] -G'oV(/? + 7V’«('*/’a,'*/’i)) + Ipi) (5.7) 


we show that H is nondecreasing in each variable. Recall Relation (3.3): 

'Ipui'lpa, V’i) = {{Si)v'<Pa - {Sa)vi’i) '>Pv{'>Pa, V’i) = {-{Si)ui’a + {Sa)u'>Pi) 
From the previous subsection, we expect |sq| and |si| to be 0{h), as well as 


\^i)u \ J \^a)u ^ ^ 

= Cl > 0 and -—— = C 2 > 0 


(5.8) 


det B det B 

Moreover, given our assumption that u = 0{pa), we have that V'0(ua) = 0, and 

Vlp{ud) = V'^(Ua) + = (5.9) 











A fast-marching algorithm for non-monotone fronts 


15 


where ^ if such that |^| = 0{h). We have: 


dH dipv (/? + 7V’«)7||^ 

(5.10) 

dlpa °yi + ( 2 ^ 7 V;„+ 72 t/; 2 +,/; 2)/^2 

(5.11) 

Wi + o{h) 

(5.12) 

= {(3-Go7)ci + 0{h) 

(5.13) 

= y 1 P Gq Cl P 0{h) > 0 for h small enough 

(5.14) 


The argument that dH/d^pi > 0 is completely analogous. We then show that the 
solver is stable in the sense of Barles & Souganidis [^. Consider again the quadratic: 

((Gq + 1)^1 ~ w‘^ 2 ((Gq + 1)A:2 ~ GqA^s) + ((^o + 1)^3 ~ Gg/^^e) = 0 

Using ( |2.6[ ), ( |3.4[ ), ( |3.8| ) and letting \si\ = Qh, where q is 0(1), we may rewrite it as: 
aw‘^ 2 {bihb2)wd{cihc2) = 0 where a, 6 i, 62 , ci, C 2 are 0(1) (5.15) 


Assuming that h < 1 , it is easily shown that the solution can be bounded by a bound 
independent of h. Verifying the consistency of the scheme reveals that it is first-order 
accurate with respect to h. Locally uniform convergence follows from [^. 

5.1.3. Iterative solver. Given a location the iterative solver returns a value 
. The initial guess is provided by the direct solver. We let: 

:= - [7 - /3^.(c,V’r)] - GA(/?+7V’.(V’^,^r))VV’2(V’^,^r) 


,,iter 


using (3.13 


where = {w'^—Wi)/\si \ for i = a^b^ c, and G'^ := F{ x(uc^, w2),t{ud, w^) )• Without 
loss of generality, we assume that the scheme is proper in the sense of 


considering 


H F ew^ if necessary. From the previous section, it follows that iL is a degenerate 
elliptic scheme. Theorem 8 of 19 guarantees convergence to the unique solution for 
arbitrary initial data. 


5.2. Global properties. 

5.2.1. The sets A and Af. We demonstrate some properties of these sets. 
Lemma 5.1. The size of JV eannot exeeed m, the number of points sampling Cq. 
Proof Initially, the size of A/" is m. Each iteration of the main loop successively: 
Removes exactly one point from JV (line 4); Adds at most one point to Af (line 12 ); 
May remove points from Af (Algorithm]^. □ 

The following lemma justifies promoting the point in Af with the smallest time 
value to the set A. 

Lemma 5.2. Subsequent iterations of the Main Loop will not affeet the value of 
the point in Af with the smallest time value. 

Proof Suppose that during the k-th iteration of the main loop, Pa £ Af is pro¬ 
moted to the set M, and let us denote by p* the point in Af with the smallest time 
value. It follows from constraint (V 2 ) that if Pa has a child point, its time value 
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satisfies td > ta At > ta > t^. Therefore pd cannot belong to the past domain of 
influence of . □ 

Remark: It is not possible to closely mimick the proofs of Sethian 


25 or Tsitsiklis 


[3^ . Indeed, those proofs assume the existence of a graph where the values of -0 (resp. 
V) are to be found. In constrast, in our framework, the choice of which point gets 
promoted to the set A directly influences the shape of the graph sampling Ai. A 

5.2.2. Global convergence. When discussing the global convergence of the 
scheme, we distinguish between different cases based on the regularity and topological 
properties of M.. 

• A4 is . Then the local solvers may be used anywhere on A4, since the 
correspondence between the LSE and the Dirichlet problem highlighted in §2.2| holds 
in the strong sense. As h ^ 0, the assumptions listed in the beginning of §5.1 
hold and imply that the local solvers are convergent. Considering the relation At = 
hjG^{Ua) + 1 as a CFL condition ensures the stability of the global scheme. Global 
convergence is thus expected. 

• A4 is and the singularities arise from topological changes. Our argument is 
that the global constraints imposed on Af are such that the local solvers are only used 
in regions of Ai that are smooth. Indeed, topological changes imply that the segments 
monitored by Book intersect. Such a situation should be prevented by Algorithm 

It is therefore expected that the scheme is convergent in those situations as well - 
although no proof is available at this point. 

• A4 is C^. Should the speed F be merely continuous, A4 may develop singulari¬ 
ties without undergoing topological changes. In those situations, it is unclear whether 


the scheme globally converges. Let us note that the iterative scheme in ^ 5.1.3 should 
converge regardless of the regularity of the Hamiltonian. Nonetheless, the size of the 
pseudo-time step At required for convergence may be hard to determine. For that 
reason, the examples considered in the next section only involve speeds that are 
everywhere along the curve. 

5.2.3. Computational time. Recall that, in what follows, n = 2 and m is 
the number of points sampling Cq. We first bound the number of operations in one 
iteration of the main loop. Finding the minimum among Af can be done in O(logm) 

Finding the L closest neighbours of Pa among 
^ mL number of operations. (Since A is ordered 


25 


iterations using a binary heap 
where L is typically ^10 requires 
in time, it is only necessary to look at the tail of the list.) 

The grid method performs at most five iterations on an grid, where s is typically 
^ 10. Computing w and evaluating constraints (VI) and (V2) at each point requires 
^ 10 operations. Evaluating constraint (E) requires ^ L operations. We arrive at a 
total of 55^^ X 10 xL = 0(10^^^) operations to find The complexity of the iterative 
solver is hard to bound a priori. Indeed, it is currently unclear how many operations 
may be required to compute the CFL condition. However, when F G x [0,T]), 

we find that devoting ^ 10^ operations to this task results in a solver that only 
requires 0(10) iterations to converge. Updating Book in Algorithmrequires 0{m) 
operations, and as noted earlier, the other procedures usually require 0(1) operations, 
but may require up to 0(m) (e.g.. When topological changes occur). We arrive at 
the conclusion that the cost of one iteration is bounded by: max{0(m), O (l0’^+^)}. 

The parameter m is determined by the user and directly influences the total 
number of points N computed by the algorithm. We have h ~ length(Co)/(2m) so 
that h = 0(I/m). 
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When F = C or F = F{t): Then M. is roughly decomposed into / strips compris¬ 
ing 0{m) points. We have I < TF/(At)min where (At)min = -1- 1 = 0{h), 

which yields a reliable estimate for the total number of points: 


A = (t^ of strips) X points per strip) 


Tf 

0{h) 


X 0{m) = 0{m^) 


(5.16) 


When F = F(x, t); Then A4 no longer exhibits a stratified structure. However, 
in general, we still expect a strip of average height At = 0{h) to contain 0{m) points, 
and therefore arrive at the same estimate for N. 

In summary, the total cost of the algorithm is GA where 

G = max{0(m), O (lO’^^^)} and A = 0{m^) (5-17) 


Remark: Although the number of operations required to find is admittedly 
large, they can be performed fairly efficiently. In Mat Lab for example, the entire 
procedure can be vectorized. Empirically, it is found to take no more than 33% of the 
total computational time. (See Appendix [D|) A 

6. Results. We illustrate the properties of the algorithm through various exam¬ 
ples. More details about the tests procedures can be found in Appendix [P} 

6.1. Examples. With the exception of Ex. (c), Co consists of a circle centered 
at the origin. The corresponding manifolds A4 are presented on Eigure | 6 Tj The exact 
normal is used to initialize the algorithm. 

faj The expanding eirele. Using the constant speed F = 1 yields a manifold Ai 
that consists of a truncated cone. 

(b) The football The speed is set to F{t) = I — resulting in a circle that 

first expands and then contracts. A4 then resembles an american football. 

(e) Two eireles. As in the previous example, the time-dependent speed F{t) = 
1 — 2t changes sign. However, the fact that Co now consists of two disjoint circles 
implies that topological changes occur during the evolution. 

(d) The oseillating eirele. The speed is set equal to F = asm{b{t -h c)) for some 
constants a, b and c. 

(e) The eseaping eirele. The speed is such that the circle first expands, and then 
moves in the positive x-direction while growing. 

(f) The 3-leaved rose. F depends on the polar angle in such a way that Ct even¬ 
tually takes the shape of a 3-leaved rose. 


6.2. Local properties. In the following, the error associated to each point is 


measured as := |0(p)| = |(/)(x, t)|, where (j) is the exact solution to the LSE (I.I). 


6.2.1. Optimization problem. We start with an analysis of the grid method 
used to solve the optimization problem. Convergence results are presented on Eigure 
6.2 for examples (a), (d) and (f). The errors associated with the solutions of the direct 
and the iterative solver are both recorded. Second order convergence with respect to 
the repulsion parameter h is clear. This result is consistent with the test procedure: 
As h gets smaller, the child point moves closer to its parents, but the distance between 
the two parent points also decreases. Eor the most general case represented by Ex. 

(f), the iterative solver slightly increases the accuracy of the solution. The number of 
iterations required to reach the tolerance is 0(1), as recorded in Table 6.1 
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(b) The football 



(c) Two circles 



(d) The oscillating circle 



(e) The escaping circle 


(f) The 3-leaved rose 



Fig. 6.1: The manifold A4 in xyt-space. points are featured in blue. 


6.2.2. Error propagation. We now discuss data obtained from running a full 
simulation. Convergence results for examples (b) and (d) are presented on Figure 
In Subfigure |6.3a| the error is seen to behave nicely even near the “tip” of A4. 


6.3 


Note that the code naturally stops once Ct = 0. Moreover, the stratified structure 
alluded to in ^5.2.3 is evident. The symmetry of Ex. (d) gives rise to cancellations 
as the circle contracts. The error appears to increase at a slow pace. 
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h (a) 

1 (d) 

(f) 

10-2 

1 

10 

9 

10-3 

1 

9 

9 

10-4 

1 

9 

8 

10-5 

4 

8 

7 

10-s 

6 

7 

6 

o 

1 

8 

10 

7 


Fig. 6.2: Local convergence results for ex- Table 6.1: Number of it- 
amples (a), (d) and (f). erations. 



(a) The football, Ex. (b) 



(b) The oscillating circle, Ex. (d)(6 periods). 


Fig. 6.3: Error versus time. 


6.3. Global properties. Define: 


L^{E) = h^Y.^p L2{E) = 




peA 


peA 


Loo{E) = maxE’p (6.1) 
peA 


with n = 2, as well as the Haussdorff distance between the reconstructed and the 
exact curve, denoted by Lh- We first consider qualitative features of the manifold Ai 
before turning to global convergence results and speed tests. 

6.3.1. Evenness of the sampling. As can be seen from Figure [6^ the sam¬ 
pling of A4 is regular. This is confirmed when the distance from a child to its parents 
is monitored in Ex. (b); see Figure 6.5 
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0 . 6 ^ 



(a) The expanding circle, Ex. (a) 



0 • pe A 

X • Desc. of P 

• Anc. of P 

★ P 

— True charac. 


(b) The expanding circle, Ex. (a) 



0.3 

0.2 

0.1 

0 

0.5 




(c) The football. Ex. (b) 


(d) The football. Ex. (b) 



V -0.5 -0.5 ^ -0.4 -0.2 0 0.2 0.4 

^ X 

(e) The 3-leaved rose. Ex. (f) (f) The 3-leaved rose. Ex. (f) 


Fig. 6.4: Sampling of A4 returned by the the algorithm, 
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(a) Distance to Parent a 


(b) Distance to Parent b 


Fig. 6.5: Three-dimensional distance from a child to each of its parents, in units of h. 


6.3.2. Domain of influence. Subfigure |6.4a may give the impression that the 
data propagate by spiralling outwards. However, this is not the case, as can be seen 
from tracking a point’s ancestry and descendance, as on Subfigures |6.4b| and |6.4d[ It 
is of particular interest to note that the true characteristic going through the point 
does lie in the numerical past and future domains of influence. 


6.3.3. Topological changes. The two-circles example illustrates the ability of 
the algorithm to deal with topological changes. As can be seen from Figure |6l6l the 
code naturally stops computing points when it reaches the ^-axis. As a result, the 
circles appropriately merge. Their separation is also well-captured. 


1.4 



(a) Full set. (b) Side view of those points with \x\ < 0.5. 


Fig. 6.6: Sampling of Ai returned by the algorithm for Ex. (c), the two circles. 

6.3.4. Convergence results. Convergence results are presented in Figure [6?f| 
along with the exact and the reconstructed curves obtained for some simulations. 
Example (e) is used to illustrate the robustness of the algorithm when M. is smooth. 
Qualitatively similar results are obtained for examples (d), (e) and (f), namely, first- 
order convergence with respect to h is observed in all the norms considered. Similar 
results are obtained for the football example, despite the fact that J\A is only C^. The 
two circles example also converges with first-order accuracy in the Li, L 2 and Ljj 
norms. Convergence in the L^c norm is not as clear. 
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(a) The escaping circle, Ex. (e). (b) The escaping circle, Ex. (e). 





(e) The two circles. Ex. (c). (up to T = 0.5) (f) The two circles. Ex. (c). 


Fig. 6. 7: (left) Global convergence results, (right) Exact and reconstructed curves. 


6.3.5. Speed tests. Our method is now compared to the standard first-order 


EMM, see for example 25 . The solution to Ex. (a) is computed on the set |x| < 0.75 


for various gridsizes - see Appendix The CPU times are presented in Subfigure 
Remark that the vertical axis is the Li norm of the error associated with the 


6.8a 


sampling of At. It is apparent that for higher accuracies, our method is faster than 
the standard EMM for this example. Subfigure |6. 8b also presents CPU times obtained 
for non-monotone examples: The trend is found to be similar to the monotone case. 
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Fig. 6.8: Accuracy of the solution vs. CPU times. 


Remark (See Appendices and [^) We point out a few changes that can be 


made to improve the computational speed of the code. As noted at the end of ^ 5.1.1 
constraints (E) and (V2) are not necessarily equivalent. However, since in practice 
(V2) (E) for an overwhelming number of points, the verification of constraint (E) 

may be taken out of the grid method, and checked a posteriori. Based on Eignre [6^ 
running the iterative solver does not necessarily improve the accuracy of the solution. 
Making the stopping criterion depend on h, e.g., — w'2\ < lOhAr may avoid 

unnecessary computations. 


Conclusion. We presented a scheme that describes the non-monotone propagation 
of fronts while featuring a numerical complexity comparable to that of the standard 
EMM. Local convergence was demonstrated and verified. Evidences of global conver¬ 
gence were supported by several examples where F G The most general case, i.e., 
F G C\ requires further theoretical investigation, and will therefore be the subject 
of subsequent papers. The theory presented in ^ and trivially extends to higher 
dimensions. Nonetheless, some practical obstacles currently prevent the design of an 
algorithm when n > 2. The most prominent one is to properly monitor the global fea¬ 
tures of the manifold, as in Algorithm]^- see umiiii for the case n = 3. Moreover, 
choosing an appropriate triplet of parents still requires some thought. 

As it stands, the algorithm is first order accurate. It may be extended to higher 
order using filtered schemes 14 32 . Allowing h to depend on (x, t) would increase 


the accuracy in regions with high mean curvature, and avoid unecessary computations 
in other regions. The precise adaptivity criteria must be carefully addressed. A 
different approach is to resort to reseeding by introducing points to the Narrow Band 
in regions of the front that are undersampled. We conclude by remarking that the 
novel ideas presented in this paper may apply to other evolution equations, such as 
linear advection, anisotropic propagation, or mean curvatuve flow. 
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Appendix A. i^^’iL’-coordinate system. 

A.l. Tilted plane. Let a unit normal vector v — ( 711 , 712 ,^ 3 ) G define a 
plane through the origin of the form t = ax^hy^ where a = —nijns and b = — 712 / 713 . 
Define y = , as well as u) = where the sign is chosen such 

that sign(u) 3 ) = sign( 773 ). The following coordinate system can be verified to satisfy 


the conditions listed in ^ 




r u = 

^ Va 2 +P' ^ ^ 

, - a , 0 ) 


V = 

\/a^+ 6 ^ ^ ^ 

, b , + 6 ^ ) 

(A.l) 

1 TC = 

dz— ( — a , 

- 6 , 1 ) 



A.2. Vertical plane. If 713 = 0, then u e describes a vertical plane. The 
orthonormal coordinate system we use is then: 


u = ( - 772 , Til , 0 ) 

u = (0,0,1) (A. 2 ) 

w = ( 77i , 772 , 0 ) 


A.3. P is the outward normal to A4. Suppose that A4 is at some point 
Pa and that = i>{pa). If \F{pa)\ > 0 , then a neighbourhood of Pa is locally described 
by '0(x), where ip satisfies [^ : 

=As 

Then 0(x, t) = sign(F(pa)) ('^(x) — t), and we have: 

u = sign(F(pa)) IL (A.4) 


Using the PDE (A.3) and the results from ^ A .2 with a = and b = we get: 


^3 = 


1 

^l + F^{Pa) 


and 


-nPa) 


(A.5) 


If F(pa) = 0, then i){pa) describes a vertical plane, and /d 3 = 1 and 73 = 0. 

Appendix B. The iterative solver. The size of the pseudo-time step used in 
the iterative solver is determined as follows. Defining the left-hand side of (3.13) as 
—H{w^)^ we set Ar = ^ ^ where 


<|(z>2-z>i).^| + |G^| 



+ |G^ -G^l 




=: Q 


e = / 7 /IO and re 1 , 7^2 G (re^ — e, re^ + e). The neighbourhood is sampled with 10 points. 

Appendix C. The grid method: Pseudo-code. This is the method we use 
in practice. Note that Lines 6-10 can be completely vectorized. 
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Algorithm 6 Constrained Optimization: The grid method 

1: h ^ h/2^ uo ^ 't^min? h + Va-i 

2 : Q = l, tol = 10 - 1 ®, i = 1 , /o = 10 ®. 

3: while Q >tol do 

4: Build a grid of 5 x 5 points, with meshsize h, centered at (i4o,r’o) 

5: On the grid, initialize / = +oo. 

6 : for each point (rt, v) on the grid do 

7: Compute w using the direct solver (using G(ua)). 

8 : Compute ^ 1 , g 2 and ^ 3 . 

9: If all three constraints are satisfied, compute f{u,v). 

10 : end for 

11 : Let the minimum of / be fi and occur at 

12: g = \fi-fi-l\ 

13: h i — ^/2, uq i — Ud-) vq i — Vd-) 'i — i T 1 

14: if i == 5 or fi == +OC then (3 = 0 end if 

15: end while 
16: Return: {ud^Vd). 


Appendix D. Examples Tests procedures. 


D.l. Examples. With the exception of the third example, Cq always consists of 
a single circle of radius vq centered at the origin. 

The expanding eirele. The speed is F = 1 and the signed distance function that 
solves the LSE is (/)(x, t) = x‘^ -\- y‘^ — t — with ro = 0.25. 

The football The speed is F{t) = l — and the signed distance function that 

solves the LSE is 0(x, t) = -s/x^F^ — (ro — (e^^ “ 1) /(ce) + 1) with ro = 0.25. 

Two eireles. The speed is F(t) = 1 — 2t and the signed distance function that 

solves the LSE up to time T = 0.5 is 0(x, t) = {x — sign(x)0.5)^ + y‘^ — (ro + t — 
with ro = 0.35. 

The oseillating eirele. The speed is F{t) = asin( 6 (t + c)) with a = 0.7, 6 = 10 
and c = 0.3, and 0(x, t) = yxFFy^ — (ro + (a/6)(cos(5c) — cos{b{t + c)))) is the 
signed distance function that solves the LSE with ro = 0.25. 

The eseaping eirele. The speed is F = {x — gt){g't-\-g)/y/{x — gty F y‘^ Fc where 
g{t) = arctan( 6 (t — 0.5)) + |, and the signed distance function that solves the LSE 
is 0(x, t) = — (ro + ct) with ro = 0.25, 6 = 10 and c = 0.5. 

The 3-leaved rose. The speed is F = cos{W) /y^l + (/t/r)^(sin(/^))^ and a solution 
of the LSE is (/)(x, y^t) = r — {t cos{W) + ro), where / = 3 is the number of petals and 
ro = 0.25. Note that this is not the signed distance function. 


D.2. Tests procedures. 

Time for eomputing The estimate of the time taken to compute provided 
in the remark at the end of ^5.2.3 is based on example (f), run with m = 150. 

Optimization problem. Each point is obtained as follows: The parent points are 
exact: Pa^ Pb ^ At. The point Pa is fixed, and pb is such that (xb^yb) = {xa^ya) + 
(—3,4)h/8. The exact normal at Pa is used to get the local coordinate frame. The 
optimization problem is solved using the direct solver and returns p^^^. Then the 


iterative solver is initialized with and run until 


w 


n+l 




< 10 -^^ Ar. 


Error propagation. Eignres [6.3a and 6.3b were produced using m = 60. 
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Global properties. The results presented on Figure [0| were obtained using M = 
25 for the expanding circle, M = 60 for the football, and M = 150 for the 3-leaved rose. 
The histograms of Figure [63] correspond to the football simulation. The binwidth is 
0.2 and the total number of points N is 6229. On Figure |6.6[ there are initially 80 
points on each circle. The data used to get the results presented on Figure [6?7| are: 
Ex. (e) Final T = 0.4, tn = 0.35, Ex. (b) tn = 0.1, Ex. (c) Einal T = 0.5, tn = 0.45. 

HausdorJJ distanee. To get the Hausdorff distance between the exact and the re¬ 
constructed curves, samplings of each set were first obtained. Eor the exact one, we 
used the exact solution to the LSE and the MatLab contour function. Eor the recon¬ 
structed one, local Delaunay triangulations were obtained using appropriate subsets 
of A. Let the resulting two cloud of points be respectively CGec and C4x- Then the 
Hausdorff distance between those sets is: 

LH(C4ec,C4x) = max < sup inf d(x,y), sup inf d(x,y)> (D.l) 

xGC^rec yGC-^ex J 

where d is the Euclidean distance function. To get the first term in braces, the exact 
signed distance function is used, since for a given x G CGec^ we have infy^^^ex <^(^5 y) = 
(I){-kAh)- To get the second one, given y G C4x, the nearest point x G C4ec is found. 

Speed tests. To get the results labelled as ‘Modified EMM’ on Eigure |6.8a[ our 
method is run. The points labelled as ‘EMM’ are obtained using the first-order East 
Marching Method with different gridsizes. Points with |x| < 0.2b-\-2dx are initialized 
with exact values. The solution to the EMM is computed on the set |x| < 0.75 + dx 
(the neighbours of those Accepted points with |x| > 0.75 were not added to the 
Narrow Band). The procedure used to sort and update the Narrow Band is the same 
in both methods, namely: • The point with the smallest time value is found using the 
Matlab min command; • This point is removed from the Narrow Band by deleting 
the corresponding row (using A7(/,:) = []); • Each new point is added at the end of 
the Narrow Band. The CPU time was evaluated using the MatLab cpu command. 
The final times of the simulations are: (a) 0.5, (c) 0.5, (d) 1.5 per., (e) 0.4, (f) 0.19. 
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